Tarea de Modelo lineales y de Sobrevivencia
Tarea 1
Parte 1
Librerias:
library(dplyr)
library(psych)
library(GGally)
library(knitr)
library(ggplot2)
library(gridExtra)
library(lmtest)
library(corrplot)
library(car)
library(datasets)
library(plotly)
library(nortest)
library(tseries)
library(carData)
library(MASS)Del siguiente link: https://rpubs.com/Joaquin_AR/226291 ,realice las siguientes acciones:
a. Replique en un Script de R, los Ejemplos 1 y 2, usando todas las funciones que se presenta el artículo.
Ejemplo 1:
Para este primer ejemplo se desea generar un modelo que permita predecir la esperanza de vida media de los habitantes de una ciudad en función de diferentes variables. Se dispone de información sobre: habitantes, analfabetismo, ingresos, esperanza de vida, asesinatos, universitarios, heladas, área y densidad poblacional. Los datos se cargan acontinuación:
datos <- as.data.frame(state.x77)
datos <- rename(habitantes = Population, analfabetismo = Illiteracy,
ingresos = Income, esp_vida = `Life Exp`, asesinatos = Murder,
universitarios = `HS Grad`, heladas = Frost, area = Area,
.data = datos)
datos <- mutate(.data = datos, densidad_pobl = habitantes * 1000 / area)1) Analizar la relación entre variables.
Acontinucación, se estudia la relación entre las variables de la base
de datos, esto para determinar si existe colinealidad o si se presentan
varibales con relaciones de tipo no lineal. Para esto se corre la
función cor() que devuelve las correlaciones entre
variables y se redondea el resultado en el tercer decimal.
## habitantes ingresos analfabetismo esp_vida asesinatos
## habitantes 1.000 0.208 0.108 -0.068 0.344
## ingresos 0.208 1.000 -0.437 0.340 -0.230
## analfabetismo 0.108 -0.437 1.000 -0.588 0.703
## esp_vida -0.068 0.340 -0.588 1.000 -0.781
## asesinatos 0.344 -0.230 0.703 -0.781 1.000
## universitarios -0.098 0.620 -0.657 0.582 -0.488
## heladas -0.332 0.226 -0.672 0.262 -0.539
## area 0.023 0.363 0.077 -0.107 0.228
## densidad_pobl 0.246 0.330 0.009 0.091 -0.185
## universitarios heladas area densidad_pobl
## habitantes -0.098 -0.332 0.023 0.246
## ingresos 0.620 0.226 0.363 0.330
## analfabetismo -0.657 -0.672 0.077 0.009
## esp_vida 0.582 0.262 -0.107 0.091
## asesinatos -0.488 -0.539 0.228 -0.185
## universitarios 1.000 0.367 0.334 -0.088
## heladas 0.367 1.000 0.059 0.002
## area 0.334 0.059 1.000 -0.341
## densidad_pobl -0.088 0.002 -0.341 1.000
Además, se representa la distribución de las vaiables mediante
histogramas. Esto con ayuda de la fucnción
multi.hist().
Otro ejemplo en el que se calculan las correlaciones y se generan
histogramas a la vez, se encuentra con la función
ggpairs().
ggpairs(datos, lower = list(continuous = "smooth"),
diag = list(continuous = "barDiag"), axisLabels = "none")De lo anterior se extra que las variables con mayor relación con
esp_vida son asesinatos,
analfabetismo y universitarios. Además, dado
que las primeras dos variables poseen cierto grado de correlación, pude
suceder que no sea útil intoducir ambas al modelo.
2) Generar el modelo.
Ahora, que se conoce la información repecto a las correlaciones y posibles distribuciones, se procede a generar un modelo mediante el método mixto y utilizando la medición AIC para la validación de los predictores.
modelo <- lm(esp_vida ~ habitantes + ingresos + analfabetismo + asesinatos +
universitarios + heladas + area + densidad_pobl, data = datos )
summary(modelo)##
## Call:
## lm(formula = esp_vida ~ habitantes + ingresos + analfabetismo +
## asesinatos + universitarios + heladas + area + densidad_pobl,
## data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.47514 -0.45887 -0.06352 0.59362 1.21823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.995e+01 1.843e+00 37.956 < 2e-16 ***
## habitantes 6.480e-05 3.001e-05 2.159 0.0367 *
## ingresos 2.701e-04 3.087e-04 0.875 0.3867
## analfabetismo 3.029e-01 4.024e-01 0.753 0.4559
## asesinatos -3.286e-01 4.941e-02 -6.652 5.12e-08 ***
## universitarios 4.291e-02 2.332e-02 1.840 0.0730 .
## heladas -4.580e-03 3.189e-03 -1.436 0.1585
## area -1.558e-06 1.914e-06 -0.814 0.4205
## densidad_pobl -1.105e-03 7.312e-04 -1.511 0.1385
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7337 on 41 degrees of freedom
## Multiple R-squared: 0.7501, Adjusted R-squared: 0.7013
## F-statistic: 15.38 on 8 and 41 DF, p-value: 3.787e-10
Con lo anterior se obtiene que el modelo pude explicar aproximadamente el 75% de la variabilidad, esto ya que el \(R^{2}\) fue de 0.7501. Por otro lado, ya que el p-valor obtenido fue significativo se pude aceptar que no se llegó al modelo por azar. Ademán, al menos uno de los coeficientes parciales de regresión es diferente de cero.
3) Selección de los mejores predictores.
Ahora se procede a seleccionar los mejores predictores empleando el método AIC.
## Start: AIC=-22.89
## esp_vida ~ habitantes + ingresos + analfabetismo + asesinatos +
## universitarios + heladas + area + densidad_pobl
##
## Df Sum of Sq RSS AIC
## - analfabetismo 1 0.3050 22.373 -24.208
## - area 1 0.3564 22.425 -24.093
## - ingresos 1 0.4120 22.480 -23.969
## <none> 22.068 -22.894
## - heladas 1 1.1102 23.178 -22.440
## - densidad_pobl 1 1.2288 23.297 -22.185
## - universitarios 1 1.8225 23.891 -20.926
## - habitantes 1 2.5095 24.578 -19.509
## - asesinatos 1 23.8173 45.886 11.707
##
## Step: AIC=-24.21
## esp_vida ~ habitantes + ingresos + asesinatos + universitarios +
## heladas + area + densidad_pobl
##
## Df Sum of Sq RSS AIC
## - area 1 0.1427 22.516 -25.890
## - ingresos 1 0.2316 22.605 -25.693
## <none> 22.373 -24.208
## - densidad_pobl 1 0.9286 23.302 -24.174
## - universitarios 1 1.5218 23.895 -22.918
## + analfabetismo 1 0.3050 22.068 -22.894
## - habitantes 1 2.2047 24.578 -21.509
## - heladas 1 3.1324 25.506 -19.656
## - asesinatos 1 26.7071 49.080 13.072
##
## Step: AIC=-25.89
## esp_vida ~ habitantes + ingresos + asesinatos + universitarios +
## heladas + densidad_pobl
##
## Df Sum of Sq RSS AIC
## - ingresos 1 0.132 22.648 -27.598
## - densidad_pobl 1 0.786 23.302 -26.174
## <none> 22.516 -25.890
## - universitarios 1 1.424 23.940 -24.824
## + area 1 0.143 22.373 -24.208
## + analfabetismo 1 0.091 22.425 -24.093
## - habitantes 1 2.332 24.848 -22.962
## - heladas 1 3.304 25.820 -21.043
## - asesinatos 1 32.779 55.295 17.033
##
## Step: AIC=-27.6
## esp_vida ~ habitantes + asesinatos + universitarios + heladas +
## densidad_pobl
##
## Df Sum of Sq RSS AIC
## - densidad_pobl 1 0.660 23.308 -28.161
## <none> 22.648 -27.598
## + ingresos 1 0.132 22.516 -25.890
## + analfabetismo 1 0.061 22.587 -25.732
## + area 1 0.043 22.605 -25.693
## - habitantes 1 2.659 25.307 -24.046
## - heladas 1 3.179 25.827 -23.030
## - universitarios 1 3.966 26.614 -21.529
## - asesinatos 1 33.626 56.274 15.910
##
## Step: AIC=-28.16
## esp_vida ~ habitantes + asesinatos + universitarios + heladas
##
## Df Sum of Sq RSS AIC
## <none> 23.308 -28.161
## + densidad_pobl 1 0.660 22.648 -27.598
## + ingresos 1 0.006 23.302 -26.174
## + analfabetismo 1 0.004 23.304 -26.170
## + area 1 0.001 23.307 -26.163
## - habitantes 1 2.064 25.372 -25.920
## - heladas 1 3.122 26.430 -23.877
## - universitarios 1 5.112 28.420 -20.246
## - asesinatos 1 34.816 58.124 15.528
##
## Call:
## lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +
## heladas, data = datos)
##
## Coefficients:
## (Intercept) habitantes asesinatos universitarios heladas
## 7.103e+01 5.014e-05 -3.001e-01 4.658e-02 -5.943e-03
Por lo tanto el mejor modelo resulta de tomar
habitantes, asesinatos,
universitarios y heladas como predictores.
modelo <- (lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +
heladas, data = datos))
summary(modelo)##
## Call:
## lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +
## heladas, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.47095 -0.53464 -0.03701 0.57621 1.50683
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.103e+01 9.529e-01 74.542 < 2e-16 ***
## habitantes 5.014e-05 2.512e-05 1.996 0.05201 .
## asesinatos -3.001e-01 3.661e-02 -8.199 1.77e-10 ***
## universitarios 4.658e-02 1.483e-02 3.142 0.00297 **
## heladas -5.943e-03 2.421e-03 -2.455 0.01802 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7197 on 45 degrees of freedom
## Multiple R-squared: 0.736, Adjusted R-squared: 0.7126
## F-statistic: 31.37 on 4 and 45 DF, p-value: 1.696e-12
Se muestra un intervalo de confianza para cada coeficiente parcial de regresión.
## 2.5 % 97.5 %
## (Intercept) 6.910798e+01 72.9462729104
## habitantes -4.543308e-07 0.0001007343
## asesinatos -3.738840e-01 -0.2264135705
## universitarios 1.671901e-02 0.0764454870
## heladas -1.081918e-02 -0.0010673977
4) Validación de condiciones para la regresión múltiple lineal.
Ahora se valida la condición de variabilidad constante para los residuos es decir la homocesdasticidad.
plot1 <- ggplot(data = datos, aes(habitantes, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot2 <- ggplot(data = datos, aes(asesinatos, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot3 <- ggplot(data = datos, aes(universitarios, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
plot4 <- ggplot(data = datos, aes(heladas, modelo$residuals)) +
geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
theme_bw()
grid.arrange(plot1, plot2, plot3, plot4)Ahora se testea la distribución normal de los residuos.
##
## Shapiro-Wilk normality test
##
## data: modelo$residuals
## W = 0.97935, p-value = 0.525
Se concluye que el análisis gráfico y la prueba, confirman la normalidad de los residuos.
Cuando graficamos los residuos en función de los valores ajustados por el modelo, esperamos ver que los residuos se distribuyan aleatoriamente alrededor de cero, sin mostrar ningún patrón. Además, esperamos que la variabilidad de los residuos sea más o menos constante a lo largo del eje X. Sin embargo, si identificamos algún patrón específico en la distribución de los residuos, como una forma cónica o una mayor dispersión en los extremos, esto indica que la variabilidad de los residuos está relacionada con los valores ajustados por el modelo, lo que sugiere una falta de homocedasticidad en el mismo.
ggplot(data = datos, aes(modelo$fitted.values, modelo$residuals)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0) +
theme_bw()##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 6.2721, df = 4, p-value = 0.1797
Por lo tanto no hay evidencia de ausencia de homocedasticidad.
Relación entre los predictores.
corrplot(cor(dplyr::select(datos, habitantes, asesinatos,universitarios,heladas)),
method = "number", tl.col = "black")Análisis de Inflación de la Varianza.
## habitantes asesinatos universitarios heladas
## 1.189835 1.727844 1.356791 1.498077
El resultado anterior señala que no hay evidencia de una correlación lineal muy alta ni inflación de varianza.
Autocorrelación.
## lag Autocorrelation D-W Statistic p-value
## 1 0.02867262 1.913997 0.818
## Alternative hypothesis: rho != 0
El resultado anterior muestra que no existe evidencia de auto correlación.
5) Identificación de posibles valores atípicos o influyentes.
Identificación de posibles valores atípicos.
datos$studentized_residual <- rstudent(modelo)
ggplot(data = datos, aes(x = predict(modelo), y = abs(studentized_residual))) +
geom_hline(yintercept = 3, color = "grey", linetype = "dashed") +
# se identifican en rojo observaciones con residuos estandarizados absolutos > 3
geom_point(aes(color = ifelse(abs(studentized_residual) > 3, 'red', 'black'))) +
scale_color_identity() +
labs(title = "Distribución de los residuos studentized",
x = "Predicción modelo") +
theme_bw() + theme(plot.title = element_text(hjust = 0.5))## integer(0)
No se encuentran valores atípicos.
## Potentially influential observations of
## lm(formula = esp_vida ~ habitantes + asesinatos + universitarios + heladas, data = datos) :
##
## dfb.1_ dfb.hbtn dfb.assn dfb.unvr dfb.hlds dffit cov.r cook.d
## Alaska 0.41 0.18 -0.40 -0.35 -0.16 -0.50 1.36_* 0.05
## California 0.04 -0.09 0.00 -0.04 0.03 -0.12 1.81_* 0.00
## Hawaii -0.03 -0.57 -0.28 0.66 -1.24_* 1.43_* 0.74 0.36
## Nevada 0.40 0.14 -0.42 -0.29 -0.28 -0.52 1.46_* 0.05
## New York 0.01 -0.06 0.00 0.00 -0.01 -0.07 1.44_* 0.00
## hat
## Alaska 0.25
## California 0.38_*
## Hawaii 0.24
## Nevada 0.29
## New York 0.23
En la tabla anterior se listan las observaciones que son significativamente influyentes en al menos uno de los predictores.
## StudRes Hat CookD
## California -0.1500614 0.38475924 0.002879053
## Hawaii 2.5430162 0.23979244 0.363778638
## Maine -2.2012995 0.06424817 0.061301962
## Nevada -0.8120831 0.28860921 0.053917754
## Washington -1.4895722 0.17168830 0.089555784
El análisis anterior muestra valores preocupantes para los valores Leverages o Distancia Cook (California y Maine). Un estudio más exhaustivo recomendaría rehacer el modelo sin dichas observaciones y comparar los resultados.
6) Conclusión.
En conclusión el modelo de regresión lineal múltiple sería: Esperanza de vida \(=5.014e^{-05}\)habitantes \(- 3.001e^{−01}\)asesinatos \(+4.658e^{−02}\)universitarios \(−5.943e^{−03}\)heladas. que es capaz de explicar el 73.6% de la variabilidad observada en la esperanza de vida.
Ejemplo 2
Este ejmplo se ejecuta mediante una base de datos con información de 30 libros. Se conoce del peso total de cada libro, el volumen que tiene y el tipo de tapas (duras o blandas). Se quiere generar un modelo lineal múltiple que permita predecir el peso de un libro en función de su volumen y del tipo de tapas.
datos <- data.frame(peso = c(800, 950, 1050, 350, 750, 600, 1075, 250, 700,
650, 975, 350, 950, 425, 725),
volumen = c(885, 1016, 1125, 239, 701, 641, 1228, 412, 953,
929, 1492, 419, 1010, 595, 1034),
tipo_tapas = c("duras", "duras", "duras", "duras", "duras",
"duras", "duras", "blandas", "blandas",
"blandas", "blandas", "blandas", "blandas",
"blandas", "blandas"))
head(datos, 4)## peso volumen tipo_tapas
## 1 800 885 duras
## 2 950 1016 duras
## 3 1050 1125 duras
## 4 350 239 duras
1) Primero se analiza las correlaciones entre la variable cualitativas.
##
## Pearson's product-moment correlation
##
## data: datos$peso and datos$volumen
## t = 7.271, df = 13, p-value = 6.262e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7090393 0.9651979
## sample estimates:
## cor
## 0.8958988
ggplot(data = datos, mapping=aes(x = tipo_tapas, y = peso, color=tipo_tapas)) +
geom_boxplot() +
geom_jitter(width = 0.1) +
theme_bw() + theme(legend.position = "none")Con los análisis anteriores se muestra que existe una relación lineal
entre la variable peso y volumen. Además la
variable tipo_tapas aparenta influir en el peso.
2) Modelo lineal multiple.
##
## Call:
## lm(formula = peso ~ volumen + tipo_tapas, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -110.10 -32.32 -16.10 28.93 210.95
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.91557 59.45408 0.234 0.818887
## volumen 0.71795 0.06153 11.669 6.6e-08 ***
## tipo_tapasduras 184.04727 40.49420 4.545 0.000672 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 78.2 on 12 degrees of freedom
## Multiple R-squared: 0.9275, Adjusted R-squared: 0.9154
## F-statistic: 76.73 on 2 and 12 DF, p-value: 1.455e-07
## 2.5 % 97.5 %
## (Intercept) -115.6237330 143.4548774
## volumen 0.5839023 0.8520052
## tipo_tapasduras 95.8179902 272.2765525
Por lo anterior, el modelo es capaz de explicar aproximadamente el 92.7% de la variabilidad observada, en el peso. Además, por el F teste, se puede decir que la muestra fue significativa.
3) Elección de los predictores.
En este caso solo hay dos posibles predicotres que mediante la
función summary() se constataron como importantes.
4) Condiciones para la regresión lineal múltiple.
Relación lineal entre los predictores numéricos y la variable dependiente:
ggplot(data = datos, aes(x = volumen, y = modelo$residuals)) +
geom_point() +
geom_smooth(color = "firebrick") +
geom_hline(yintercept = 0) +
theme_bw()Se satisface la condición de linealidad y se observa un dato que podría ser atípico.
Distribución normal de los residuos:
##
## Shapiro-Wilk normality test
##
## data: modelo$residuals
## W = 0.85497, p-value = 0.02043
Por lo anterior, se concluye que la condición de normalidad no se cumple, posiblemente por el dato atípico observado anteriormente. Por lo tanto, se vuelve a correr el código sin el valor atípico.
##
## Shapiro-Wilk normality test
##
## data: modelo$residuals[-dato_at]
## W = 0.9602, p-value = 0.7263
Por lo tanto, los residuos sí se distribuyen de forma normal.
Variabilidad constante de los residuos:
ggplot(data = data.frame(predict_values = predict(modelo),
residuos = residuals(modelo)),
aes(x = predict_values, y = residuos)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0) +
theme_bw()##
## studentized Breusch-Pagan test
##
## data: modelo
## BP = 2.0962, df = 2, p-value = 0.3506
No hay evidencias de falta de homocedasticidad. Además, dado que solo hay un predictor cuantitativo, no se puede dar colinealidad.
Autocorrelación:
## lag Autocorrelation D-W Statistic p-value
## 1 0.0004221711 1.970663 0.744
## Alternative hypothesis: rho != 0
No hay evidencias de autocorrelación.
5) Identificación de posibles valores atípicos o influyentes
## rstudent unadjusted p-value Bonferroni p
## 13 5.126833 0.00032993 0.004949
Justo como se vio anteriormente, existe un dato atípico.
## Potentially influential observations of
## lm(formula = peso ~ volumen + tipo_tapas, data = datos) :
##
## dfb.1_ dfb.vlmn dfb.tp_t dffit cov.r cook.d hat
## 4 -0.16 0.18 -0.10 -0.23 1.98_* 0.02 0.36
## 11 0.70 -1.26_* 0.57 -1.54_* 0.83 0.64 0.38
## 13 0.31 0.67 -1.31_* 2.07_* 0.04_* 0.46 0.14
## StudRes Hat CookD
## 4 -0.3008849 0.3616894 0.01850173
## 11 -1.9897106 0.3757842 0.63729788
## 13 5.1268330 0.1397761 0.45819722
El análisis anterior muestra la existencia de múltiples observaciones influyentes, pero ninguna preocupante bajo los valores leverage hat o distancia Cook.
6) Conclusión.
El modelo Peso \(=13.91557+0.71795\)volumen \(+184.04727\)tipotapas, es capaz de explicar el 92.7% de la variabilidad observada en el pseo de los libros.
b.Explique brevemente las siguientes funciones, incluyendo valores recibe y cuáles son sus salidas.
multi.histde la bibliotecapsych: Recibe una matriz o dataframe y devuelve un plot con multiples histogramas para cada variable. Además, puede representar distribuciones de densidad para cada gráfico.ggpairsde la bibliotecaGGally: Esta función recibe un dataframe con variables continuas o categoricas y devuleve un gráfico que contiene las correlaciones encima de la diagonal, por debajo los diagramas de dispersión entre las variables continuas, en la diagonal los gráficos de densidad de las variables continuas y en los lados los histogramas y box plots de la combinación entre variables continuas y categóricas.- Función
Stepde la sección 3 del ejemplo 1: Esta función permite encontrar el mejor modelo basado en AIC utilizando a elección el método forward, backward o moxto. Esta función recibe un modelo el cual es un objeto de tipolmoglm, y devuelve el modelo seleccionado. Confint: Recibe un objeto de tipo modelo y devuelve una matriz o vector con los limites superior e inferior de confianza para cada parámetro.Corrplotde la bibliotecacorrplot: Esta función recibe una matriz de correlaciones y devuelve un gráfico de dichas correlaciones.Dwtde la bibliotecacar: Esta función recibe una serie de tiempo univariada o multivariada (vectores numéricos, matrices y dataframes son también aceptados) y calcula los coeficientes de la transformada wavelet discretaDevuelve. Devuelve un objeto de tipodwtcon información respecto a la autocorrelación.InfluencePlot: Esta función recibe un modelo lineal y retorna un gráfico de “burbujas” de los residuos de Studentized frente a los valores hat, con las áreas de los círculos que representan las observaciones proporcionales al valor Distancia de Cook.outlierTestde la bibliotecacar: Esta función recibe un modelo de regresión lineal e informa de los valores atípicos basado en los p-valores de Bonferroni.influence.measure: Esta función se utiliza para calcular algunos de los diagnósticos de regresión para modelos lineales y lineales generalizados. Recibe un modelo lineal y retorna objeto que al pasarlo por la funciónsummary, revela información relevante sobre la influencia de las observaciones.
c.Explique los test incluyendo cual es la hipótesis nula.
Shapiro-Wilk para normalidad: Esta prueba estadística se utiliza para determinar si un conjunto de datos proviene de una distribución normal, para ello se establecen las hipótesis: \(H_{0}:\) La distribución es normal y \(H_{1}:\) La distribución no es normal. Al realizar la prueba y dado un nivel de significancia, el p-valor arrojado indicará si se rechaza la hipotesis nula (\(H_{0}\)) o sí, por lo contrario, se acepta.
Studentized Breusch-Pagan para homocedasticidad: Esta prueba se implementa para corroborar la homocedasticidad al analizar si la varianza estimada de los residuos de una regresión dependen de los valores de las variables independientes. En este test las hipótesis son: \(H_{0}:\) Los errores tiene varianza constante y \(H_{1}:\) Los errores no tiene varianza constante. Además, como en el test anterior dado un nivel de significancia, el p-valor arrojado podrá o no rechazar la hipótesis nula.
Parte 2
## Population Income Illiteracy Life Exp Murder HS Grad Frost Area
## Alabama 3615 3624 2.1 69.05 15.1 41.3 20 50708
## Alaska 365 6315 1.5 69.31 11.3 66.7 152 566432
## Arizona 2212 4530 1.8 70.55 7.8 58.1 15 113417
## Arkansas 2110 3378 1.9 70.66 10.1 39.9 65 51945
## California 21198 5114 1.1 71.71 10.3 62.6 20 156361
## Colorado 2541 4884 0.7 72.06 6.8 63.9 166 103766
1.Cálcule la matriz de Correlaciones de la base de datos states.x77
## Population Income Illiteracy Life Exp Murder
## Population 1.00000000 0.2082276 0.10762237 -0.06805195 0.3436428
## Income 0.20822756 1.0000000 -0.43707519 0.34025534 -0.2300776
## Illiteracy 0.10762237 -0.4370752 1.00000000 -0.58847793 0.7029752
## Life Exp -0.06805195 0.3402553 -0.58847793 1.00000000 -0.7808458
## Murder 0.34364275 -0.2300776 0.70297520 -0.78084575 1.0000000
## HS Grad -0.09848975 0.6199323 -0.65718861 0.58221620 -0.4879710
## Frost -0.33215245 0.2262822 -0.67194697 0.26206801 -0.5388834
## Area 0.02254384 0.3633154 0.07726113 -0.10733194 0.2283902
## HS Grad Frost Area
## Population -0.09848975 -0.3321525 0.02254384
## Income 0.61993232 0.2262822 0.36331544
## Illiteracy -0.65718861 -0.6719470 0.07726113
## Life Exp 0.58221620 0.2620680 -0.10733194
## Murder -0.48797102 -0.5388834 0.22839021
## HS Grad 1.00000000 0.3667797 0.33354187
## Frost 0.36677970 1.0000000 0.05922910
## Area 0.33354187 0.0592291 1.00000000
2. Genere tres modelos lineales usando la función lm de R para:
Modelo 1: Murder ~ Population
##
## Call:
## lm(formula = Murder ~ Population, data = state.x77)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9855 -3.0119 -0.3128 2.4986 7.9014
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.1713934 0.6869410 8.984 7.49e-12 ***
## Population 0.0002841 0.0001121 2.535 0.0146 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.503 on 48 degrees of freedom
## Multiple R-squared: 0.1181, Adjusted R-squared: 0.09972
## F-statistic: 6.427 on 1 and 48 DF, p-value: 0.01455
Modelo 2: Murder ~ Illiteracy
##
## Call:
## lm(formula = Murder ~ Illiteracy, data = state.x77)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5315 -2.0602 -0.2503 1.6916 6.9745
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.3968 0.8184 2.928 0.0052 **
## Illiteracy 4.2575 0.6217 6.848 1.26e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.653 on 48 degrees of freedom
## Multiple R-squared: 0.4942, Adjusted R-squared: 0.4836
## F-statistic: 46.89 on 1 and 48 DF, p-value: 1.258e-08
Modelo 3: Murder ~ Population + Illiteracy
##
## Call:
## lm(formula = Murder ~ Population + Illiteracy, data = state.x77)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7652 -1.6561 -0.0898 1.4570 7.6758
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.652e+00 8.101e-01 2.039 0.04713 *
## Population 2.242e-04 7.984e-05 2.808 0.00724 **
## Illiteracy 4.081e+00 5.848e-01 6.978 8.83e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.481 on 47 degrees of freedom
## Multiple R-squared: 0.5668, Adjusted R-squared: 0.5484
## F-statistic: 30.75 on 2 and 47 DF, p-value: 2.893e-09
3. Represente gráficamente los gráficos de dispersión correspondientes, las rectas de mínimos cuadrados para los tres modelos anteriores
de mínimos cuadrados para los tres modelos anteriores.
# Gráfico de dispersión y recta de mínimos cuadrados para el modelo 1: Murder ~ Population
plot(state.x77$Population, state.x77$Murder, main = "Murder vs Population", xlab = "Population", ylab = "Murder")
abline(model1, col = "red")# Gráfico de dispersión y recta de mínimos cuadrados para el modelo 2: Murder ~ Illiteracy
plot(state.x77$Illiteracy, state.x77$Murder, main = "Murder vs Illiteracy", xlab = "Illiteracy", ylab = "Murder")
abline(model2, col = "blue")# Gráfico de dispersión y recta de mínimos cuadrados para el modelo 3: Murder ~ Population + Illiteracy
# Crear una cuadrícula de valores para Illiteracy y Population
illiteracy_grid <- seq(min(state.x77$Illiteracy), max(state.x77$Illiteracy), length.out = 20)
population_grid <- seq(min(state.x77$Population), max(state.x77$Population), length.out = 20)
grid <- expand.grid(Illiteracy = illiteracy_grid, Population = population_grid)
# Calcular los valores predichos de Murder para cada combinación de Illiteracy y Population
predictions <- predict(model3, newdata = grid)
# Calcular los valores predichos de Murder para los datos originales
predictions_original <- predict(model3)
# Coeficientes del plano de mínimos cuadrados
intercept <- coef(model3)[1]
slope_population <- coef(model3)[2]
slope_illiteracy <- coef(model3)[3]
# Calcular los valores de Murder para la recta de mínimos cuadrados
murder_line <- intercept + slope_population * grid$Population + slope_illiteracy * grid$Illiteracy
# Crear el gráfico 3D con plotly
scatter_plot <- plot_ly(x = grid$Illiteracy, y = grid$Population, z = predictions, type = "surface")
scatter_plot <- add_trace(scatter_plot, x = state.x77$Illiteracy, y = state.x77$Population, z = state.x77$Murder, mode = "markers", type = "scatter3d", marker = list(size = 5))
# Agregar plano de mínimos cuadrados
scatter_plot <- add_trace(scatter_plot, x = grid$Illiteracy, y = grid$Population, z = murder_line, mode = "lines", type = "scatter3d", line = list(color = "blue", width = 3))
# Etiquetas y título
scatter_plot <- layout(scatter_plot, scene = list(xaxis = list(title = "Illiteracy"), yaxis = list(title = "Population"), zaxis = list(title = "Murder")), title = "Murder vs Population + Illiteracy")
# Mostrar el gráfico
scatter_plot4. Evalue para el modelo Murder ~ Population las hipótesis de homocedasticidad, normalidad de errores y ausencia de puntos influyentes o atípicos.
Homocedasticidad
##
## studentized Breusch-Pagan test
##
## data: model1
## BP = 1.4303, df = 1, p-value = 0.2317
Como el valor p (0.2317) es mayor que el nivel de significancia típico de 0.05, no hay suficiente evidencia para rechazar la hipótesis nula de homocedasticidad.Por ende, según esta prueba, el modelo Murder ~ Population cumple con el supuesto de homocedasticidad.
Normalidad de errores
##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.95662, p-value = 0.0642
Dado que el valor p (0.0642) es mayor que el nivel de significancia típico de 0.05, no hay suficiente evidencia para rechazar la hipótesis nula de normalidad. Por lo tanto, según esta prueba, los residuos del modelo parecen seguir una distribución normal.
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuos
## D = 0.080377, p-value = 0.5801
Dado que el valor p (0.5801) es mayor que el nivel de significancia comúnmente utilizado de 0.05, no hay suficiente evidencia para rechazar la hipótesis nula de normalidad. Por lo tanto, según esta prueba, los residuos del modelo parecen seguir una distribución normal.
##
## Jarque Bera Test
##
## data: residuos
## X-squared = 2.603, df = 2, p-value = 0.2721
Dado que el valor p (0.2721) es mayor que el nivel de significancia típicamente utilizado de 0.05, no hay suficiente evidencia para rechazar la hipótesis nula de que los datos provienen de una distribución normal. Por lo tanto, en este caso, los datos podrían considerarse aproximadamente normales según el test de Jarque-Bera.
Puntos influyentes o atípicos
num_bins <- 30
hist(residuos, breaks = num_bins, col = "skyblue", main = "Histograma de Residuos", xlab = "Residuos")mediaRes <- mean(residuos)
desRes <- sd(residuos)
valorRefAtipD <- mediaRes+2*desRes
valorRefAtipI <- mediaRes-2*desRes
atipDerecha <- residuos[residuos > valorRefAtipD]
atipIzquierda <- residuos[residuos < valorRefAtipI]
# Imprime los residuos atípicos
print("Residuos atípicos en el extremo derecho:")## [1] "Residuos atípicos en el extremo derecho:"
## Alabama
## 7.901416
## [1] "Residuos atípicos en el extremo izquierdo:"
## named numeric(0)
# Prueba de los residuos estandarizados
residuos_estandarizados <- residuos / sd(residuos)
valores_atipicos_residuos <- which(abs(residuos_estandarizados) > 2) # Umbral común: 2 desviaciones estándar
print(valores_atipicos_residuos)## Alabama
## 1
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## Alabama 2.388285 0.020997 NA
Según los gráficos y pruebas anteriores, hay evidencia de al menos un
valor atípico (Alabama) por la derecha.
## StudRes Hat CookD
## Alabama 2.3882848 0.02040822 0.054112891
## California -0.6492747 0.31422549 0.097758047
## Georgia 1.8723014 0.02047985 0.034828716
## New York -0.1300681 0.21582941 0.002376835
Además según lo anteior, se presentan ciertos valores influyentes.
5. Estima las tasas de homicidio para niveles de Illiteracy de 0.25, 1.2, 2.1, 3.0, 4.0. en el segundo modelo.
coeficientes <- coef(model2)
niveles_illiteracy <- c(0.25, 1.2, 2.1, 3.0, 4.0)
tasas_homicidioE <- coeficientes[1] + coeficientes[2] * niveles_illiteracy
resultado <- data.frame(Illiteracy = niveles_illiteracy, Tasa_Homicidio_Estimada = tasas_homicidioE)
print(resultado )## Illiteracy Tasa_Homicidio_Estimada
## 1 0.25 3.461140
## 2 1.20 7.505724
## 3 2.10 11.337435
## 4 3.00 15.169146
## 5 4.00 19.426603
Parte 3
Se carga la base de datos correspondiente
Y = c(47.83,59.51,50.38,53.24,47.26,50.52)
Ynames = c("obs1","obs2","obs3","obs4","obs5","obs6")
names = c("b0", "tasaFlujoGas1", "tasaFlujoGas2", "aperturaBoq", "tempGas")
X = matrix(c(1, 124.17,134.07,23.32,210.11,
1,149.46,142.21,41.82,229.67,
1,131.65,146.62,21.14,231.10,
1,139.49,136.16,45.79,206.03,
1,113.03,125.41,41.51,222.67,
1,134.57,165.84,32.42,219.59),nrow=6,byrow = TRUE,dimnames=list(Ynames,names))
X2 = data.frame(cbind(Y,X[,-1]))a. Estime un modelo de regresión lineal, con la función lm() para la variable X2.
##
## Call:
## lm(formula = Y ~ ., data = X2)
##
## Residuals:
## obs1 obs2 obs3 obs4 obs5 obs6
## 0.53326 0.48658 -0.67440 -0.57468 0.03905 0.19019
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.23386 12.59988 -0.812 0.5657
## tasaFlujoGas1 0.33540 0.05226 6.417 0.0984 .
## tasaFlujoGas2 -0.07329 0.04757 -1.541 0.3665
## aperturaBoq 0.08865 0.05935 1.494 0.3756
## tempGas 0.11253 0.05326 2.113 0.2814
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.159 on 1 degrees of freedom
## Multiple R-squared: 0.9867, Adjusted R-squared: 0.9334
## F-statistic: 18.51 on 4 and 1 DF, p-value: 0.1724
b. La función lm() estima la inversa de XtX con base en el proceso de descomposición QR, de Gram–Schmidt, realice un breve descripción de este proceso y verifique usando la función qr.solve () que los estimadores del resultado de la regresión lineal de punto a son los mismos.
El proceso de descomposición QR se desarrolla de la siguiente forma:
Sea \(X\) una matriz nxm . Si consideramos las m columnas de \(X\), denotadas por \(x_j\) podemos obtener una base ortogonal de vectores unitarios \(q_j\) tales que:
Si \(X\) tiene rango completo (es decir, sus columnas son linealmente independientes), podemos relacionar los vectores \(x_j\) y \(q_j\) mediante la ecuación
\(x_k = r_{1k}q_1 + r_{2k}q_2 + r_{3k}q_3 + ... + r_{kk}q_k.\), con k = 1, …,m
Como los vectores \(q_j\) son ortogonales, los coeficientes vienen dados por
\(rij = q_i*a_j\) con \(i \not = j\)
Esto permite escribir
\(X = QR\) , donde \(Q\) es una matriz nxm que satisface que \(QQ = I\). Y, \(R\) es una matriz mxm triangular superior.
Si \(X\) no tiene rango completo, entonces basta con seleccionar \(q_k\) de manera tal que sea ortogonal a \((q_1, . . . , q_{k-1})\). De esta manera también se puede escribir \(X = QR\).
Finalmente, la inversa de la matriz \(X^tX\) se puede calcular mediante la descomposición \(QR\) de la matriz \(X\).
Verificación usando qr.solve()
Empleando la función qr.solve se obtiene lo siguiente:
## b0 tasaFlujoGas1 tasaFlujoGas2 aperturaBoq tempGas
## -10.23386366 0.33540014 -0.07329371 0.08864733 0.11252772
Como se puede observar, son los mismos valores de beta que brinda la función lm en el punto a.
c. Usando una matriz pseudo-inversa (Moonre-Penrose) y según lo visto en clases, reconstruya los valores de la regresión del punto a.
Estimadores
Primeramente, se calcula la matriz ortogonal
## [,1] [,2] [,3] [,4] [,5]
## [1,] 118.13156142 -0.0340693673 -0.1000545541 -0.1742167735 -0.4243787993
## [2,] -0.03406937 0.0020325502 -0.0010235978 -0.0010630811 -0.0002401096
## [3,] -0.10005455 -0.0010235978 0.0016840371 0.0009228886 -0.0001597012
## [4,] -0.17421677 -0.0010630811 0.0009228886 0.0026209996 0.0004267736
## [5,] -0.42437880 -0.0002401096 -0.0001597012 0.0004267736 0.0021107277
La solución de los betas es la siguiente:
## [,1]
## [1,] -10.23386366
## [2,] 0.33540014
## [3,] -0.07329371
## [4,] 0.08864733
## [5,] 0.11252772
Los cuales son los mismo que los obtenidos mediante la función lm.
Std error de los estimadores
Como primer paso, se calcula la matriz de covarianzas
## (Intercept) tasaFlujoGas1 tasaFlujoGas2 aperturaBoq
## (Intercept) 158.75700508 -0.0457858226 -0.1344633150 -0.234129921
## tasaFlujoGas1 -0.04578582 0.0027315443 -0.0013756131 -0.001428675
## tasaFlujoGas2 -0.13446331 -0.0013756131 0.0022631774 0.001240270
## aperturaBoq -0.23412992 -0.0014286747 0.0012402700 0.003522361
## tempGas -0.57032267 -0.0003226833 -0.0002146225 0.000573541
## tempGas
## (Intercept) -0.5703226673
## tasaFlujoGas1 -0.0003226833
## tasaFlujoGas2 -0.0002146225
## aperturaBoq 0.0005735410
## tempGas 0.0028366070
Posteriormente, se obtiene las std error de los betas
errores = c()
for (i in 1: 5) {
errores[i] = sqrt(Var_Betas[i,i])
}
errores_Betas = matrix(errores, dimnames = list(names, "Std error" ))
errores_Betas## Std error
## b0 12.59988115
## tasaFlujoGas1 0.05226418
## tasaFlujoGas2 0.04757286
## aperturaBoq 0.05934948
## tempGas 0.05325981
Como se puede observar, son iguales a los errores calculados con la
función lm.
T-value
Se procede a calcular el los T-valores de cada variable:
t_valor_matriz = c()
for (i in 1: 5) {
t_valor_matriz[i] = Betas[i,]/errores_Betas[i,]
}
t_valores = matrix(t_valor_matriz,dimnames = list(names, "T-value" ))
t_valores## T-value
## b0 -0.8122191
## tasaFlujoGas1 6.4174002
## tasaFlujoGas2 -1.5406622
## aperturaBoq 1.4936495
## tempGas 2.1128076
Pr(>|t|)
Con lo anterior, se obtiene que las probabilidades son:
probabilidades = c() for (i in 1: 5) { probabilidades[i] = 2*(1-pt(abs(t_valor_matriz[i]), 1)) } P_T = matrix(probabilidades,dimnames = list(names, "Pr(>|t|)" )) P_T## Pr(>|t|) ## b0 0.56573154 ## tasaFlujoGas1 0.09841069 ## tasaFlujoGas2 0.36651617 ## aperturaBoq 0.37558169 ## tempGas 0.28142638De tal manera, comparando con lo obtenido por lm, se tienen los mismos resultados.
Residual estándar error.
Primero, se obtiene la matriz de proyección P
## obs1 obs2 obs3 obs4 obs5 obs6
## obs1 0.78840211 -0.19307656 0.26760326 0.22803272 -0.015495472 -0.075466065
## obs2 -0.19307656 0.82382358 0.24417974 0.20807284 -0.014139141 -0.068860462
## obs3 0.26760326 0.24417974 0.66156796 -0.28838804 0.019596787 0.095440297
## obs4 0.22803272 0.20807284 -0.28838804 0.75425595 0.016699007 0.081327525
## obs5 -0.01549547 -0.01413914 0.01959679 0.01669901 0.998865255 -0.005526437
## obs6 -0.07546607 -0.06886046 0.09544030 0.08132752 -0.005526437 0.973085142
Luego, se calcula la matriz IP
Tercero, se obtiene la SSE (suma de los cuadrados residuales)
Finalmente, se calcula el error residual estándar
## [,1]
## [1,] 1.159267
La cual es la misma que da la función lm.
R cuadrado
Primeramente, se obtiene la SST (total de las sumas la cuadrado)
## [,1]
## [1,] 15987.57
Se centran las observaciones, por lo que se tiene el SSTC (SST
centrado)
## [,1]
## [1,] 100.8377
Por consiguiente, el R cuadrado es:
## [,1]
## [1,] 0.9866726
que es el mismo que el dado por la función lm.
- R cuadrado Ajustado
## [,1]
## [1,] 0.9333632
Es posible identificar que es igual al calculado por la función lm.
F-estadistico y su p-value.
Para obtener el F-estadístico, es necesario conocer sus grados de libertad. Se sabe que posee t y N-r = 1 grados de libertad. Entonces, queda pendiente conocer t.
Por medio de las probabilidades Pr(>|t|), se puede identificar que la variable tasaFlujoGas1 posee la menor probabilidad y por ende se considera como una variable significativa junto a B0, con respecto a las otras. Dado que t es la cantidad de variables que no aportan información significativa, se tiene que t = 4.
De tal manera, el estadístico F es:
## [,1]
## [1,] 18.50841
Finalmente, el p-valor corresponde a :
## [,1]
## [1,] 0.1723969
Como se puede ver,ambos valores son los mismos que los dados por
la función lm.
d. Realice una prueba de contraste para determinar si se acepta o no la hipótesis que Tasa Flujo Gas 2 + aperturaBoq = 0.
Se realiza una prueba de contraste con las siguientes hipótesis:
- \(H_0 : B2 + B3 = 0\)
- \(H_1 : B2 + B3 \not = 0\)
Primero, se calcular el error:
## [1] 0.09091798
Así, el T-valor es:
## [1] 0.1689435
Finalmente, se calcula el p-valor para determinar si se rechaza o no la hipótesis nula. Para este, tenemos r = 5 y N = 6. De tal manera, el p-valor es:
## [1] 0.8934533
Por tanto, con un nivel de confianza del 95%, se acepta la hipótesis nula.